* ********************************************************************************************
* The Effects of Weather Shocks on Economic Activity: What are the Channels of Impact?
* Sebastian Acevedo, Mico Mrkaic, Natalija Novta, Evgenia Pugacheva, Petia Topalova
* With support from Gavin Asdorian, Olivia Ma, Jilun Xing and Yuan Zeng 
* Replication files for Table 1, column 9
* ********************************************************************************************
clear all
set more off
set matsize 11000

cd "C:\Users\..\replication" // set your working directory to where the replication folder is saved

* the Y variable that will be used in regressions: ngdprpc
global Y ngdprpc

* Horizons
global k 7

* ****************************************************************************************
* Data
* ****************************************************************************************
use "data/subnational_dataset.dta", clear

drop if year > 2016

foreach X in pwtemp pwprecip {
    gen `X'2=`X'^2
}

gen sqyear=year^2

by provcode (year), sort: gen d0_ln_gdppc = ln_gdppc - l.ln_gdppc
forval t = 1/$k {
    by provcode (year), sort: gen d`t'_ln_gdppc = f`t'.ln_gdppc - l.ln_gdppc
}

* Create WDI region year fixed effects
gen tempvar1 = wdi_region + " " + string(year) if wdi_region != ""
encode tempvar1, generate(ry)
drop tempvar1

xi i.ifscode*year i.ifscode*sqyear  i.ry i.year


* ****************************************************************************************
* Variables where to store results
* ****************************************************************************************

gen period = -1 in 1

foreach X in t t2 p p2 t_ae t_em t_lic p_ae p_em p_lic{
	gen `X'_coef = 0 in 1
	gen `X'_se = 0 in 1
	gen `X'_pvalue = 0 in 1
}


gen nobs = 0 in 1
gen ncountry = 0 in 1
gen nprov = 0 in 1
gen r2 = 0 in 1

* ****************************************************************************************
* Program
* ****************************************************************************************

forval t = 0/$k {
	areg d`t'_ln_gdppc pwtemp pwtemp2 pwprecip pwprecip2 l.pwtemp l.pwtemp2 l.pwprecip l.pwprecip2 l.d0_ln_gdppc _Iry*, cluster(provcode) absorb(provcode)
	
	replace period = `t' in `=`t'+2'

	* temperature
	replace t_coef = _b[pwtemp] in `=`t'+2'
	replace t_se = _se[pwtemp] in `=`t'+2'
	replace t_pvalue = 2 * ttail(e(df_r),abs(_b[pwtemp]/_se[pwtemp])) in `=`t'+2'

	* temperature squared
	replace t2_coef = _b[pwtemp2] in `=`t'+2'
	replace t2_se = _se[pwtemp2] in `=`t'+2'
	replace t2_pvalue = 2 * ttail(e(df_r),abs(_b[pwtemp2]/_se[pwtemp2])) in `=`t'+2'

	* precipitation
	replace p_coef = _b[pwprecip] in `=`t'+2'
	replace p_se = _se[pwprecip] in `=`t'+2'
	replace p_pvalue = 2 * ttail(e(df_r),abs(_b[pwprecip]/_se[pwprecip])) in `=`t'+2'

	* precipitation squared
	replace p2_coef = _b[pwprecip2] in `=`t'+2'
	replace p2_se = _se[pwprecip2] in `=`t'+2'
	replace p2_pvalue = 2 * ttail(e(df_r),abs(_b[pwprecip2]/_se[pwprecip2])) in `=`t'+2'

	replace nobs = `e(N)' in `=`t'+2'
	replace r2 = `e(r2_a)' in `=`t'+2'

	* number of provinces
	replace nprov = `e(N_clust)' in `=`t'+2'

	* number of countries
	levelsof ifscode if e(sample), local(country_list)
	local country_count: list sizeof country_list
	replace ncountry = `country_count' in `=`t'+2'

	*** Temperature
	* lincom AE
	lincom _b[pwtemp] + (2 * _b[pwtemp2] * 11)
	replace t_ae_coef = r(estimate) in `=`t'+2'
	replace t_ae_se = r(se) in `=`t'+2'
	replace t_ae_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'

	* lincom EM
	lincom _b[pwtemp] + (2 * _b[pwtemp2] * 22)
	replace t_em_coef = r(estimate) in `=`t'+2'
	replace t_em_se = r(se) in `=`t'+2'
	replace t_em_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'

	* lincom LIC
	lincom _b[pwtemp] + (2 * _b[pwtemp2] * 25)
	replace t_lic_coef = r(estimate) in `=`t'+2'
	replace t_lic_se = r(se) in `=`t'+2'
	replace t_lic_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'

	*** Precipitation
	* lincom AE
	lincom _b[pwprecip] + (2 * _b[pwprecip2] * 8)
	replace p_ae_coef = r(estimate) in `=`t'+2'
	replace p_ae_se = r(se) in `=`t'+2'
	replace p_ae_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'

	* lincom EM
	lincom _b[pwprecip] + (2 * _b[pwprecip2] * 9)
	replace p_em_coef = r(estimate) in `=`t'+2'
	replace p_em_se = r(se) in `=`t'+2'
	replace p_em_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'

	* lincom LIC
	lincom _b[pwprecip] + (2 * _b[pwprecip2] * 11)
	replace p_lic_coef = r(estimate) in `=`t'+2'
	replace p_lic_se = r(se) in `=`t'+2'
	replace p_lic_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
}

export excel period - r2 in 1/9 using "output/Table_1.xlsx", sheet("column_9") sheetreplace firstrow(var)

